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We describe a first-principles method for calculating electronic structure, vibrational modes and 
frequencies, electron-phonon couplings, and inelastic electron transport properties of an atomic- 
scale device bridging two metallic contacts under nonequilibrium conditions. The method extends 
the density-functional codes Siesta and Transiesta that use atomic basis sets. The inelastic 
conductance characteristics are calculated using the nonequilibrium Green's function formalism, 
and the electron-phonon interaction is addressed with perturbation theory up to the level of the 
self-consistent Born approximation. While these calculations often are computationally demanding, 
we show how they can be approximated by a simple and efficient lowest order expansion. Our method 
also addresses effects of energy dissipation and local heating of the junction via detailed calculations 
of the power flow. We demonstrate the developed procedures by considering inelastic transport 
through atomic gold wires of various lengths, thereby extending the results presented in [Frederiksen 
et al. Phys. Rev. Lett. 93, 256601 (2004)]. To illustrate that the method applies more generally to 
molecular devices, we also calculate the inelastic current through different hydrocarbon molecules 
between gold electrodes. Both for the wires and the molecules our theory is in quantitative agreement 
with experiments, and characterizes the system-specific mode selectivity and local heating. 

PACS numbers: 63.22.-|-m, 71.15.-m, 72.10.-d, 73.23.-b 



I. INTRODUCTION 

Electron transport in atomic-scale devices is an 
important research area where both fundamental 
physics and technological opportunities are simultane- 
ously addressed ii Examples of novel structures include 
molecules in self-assembled monolayers (SAM)^ carbon 
nanotube based components?^ nanowireSf^ and single- 
molecule junctions.— '^li'^i'^ Also conventional lithography- 
based semiconductor electronics is rapidly being pushed 
towards the scale where atomic features become impor- 
tant. For example, the transistor gate oxide is now only 
a few atomic layers thickJ^ 

The interaction between electrons and nuclear vibra- 
tions plays an important role for the electron transport at 
the nanometer scalcj ^^d^ and is being addressed experi- 
mentally in ultimate atomic-sized systemsi -'^^d^'-'^^d^d'^d^ 
Effects on the electronic current due to energy dissipa- 
tion from electron-phonon (e-ph) interactions are rele- 
vant, not only because they affect device characteristics, 
induce chemical reactions)^ and ultimately control the 
stability; these may also be used for spectroscopy to de- 
duce structural information — such as the bonding con- 
figuration in a nanoscale junction — which is typically not 
accessible by other techniques simultaneously with trans- 
port measurements. 

The signatures of e-ph interaction have been observed 
in a variety of nanosystems. In the late 1990s in- 
elastic electron tunneling spectroscopy (lETS) on single 
molecules was successfully demonstrated using a scan- 
ning tunneling microscope (STM).— Later, in the quan- 
tum dot regime, measurements on a single Cgo transis- 
tor showed features indicating a strong coupling between 



center-of-mass motion of the molecule and single-electron 
hoppingii^ Point contact spectroscopy has also revealed 
phonon signals in the high-conductance regime, e.g., in 
atomic wiresi^ and individual molecules Most recently, 
inelastic measurements have also been reported on SAMs 
of alkyl- and 7r-conjugated molecular wiresi^iSii^ These 
developments show the need for fully atomistic quantita- 
tive theories to accurately model structural, vibrational, 
and transport properties of nanoscale systems. 

The density functional theory (DFT) approach offers 
an atomistic description of total energy properties of 
nanosystems without system specific adjustable param- 
eters. Furthermore, in combination with the nonequi- 
librium Green's function (NEGF) methodS^i^ it has re- 
cently become a popular approach to quantum trans- 
port in atomic structuresj ^^'^^i^''''^^i^^'^°i'^-^i'^^ From the 
comparison with experimental data it has been estab- 
lished that total energy properties such as atomic struc- 
ture and vibrations in general arc well described by 
DFT with the local or gradient approximations for ex- 
change and correlationi^ However, while transport prop- 
erties may also be calculated from DFT this is not rig- 
orously justifiedi^i^ On the other hand such an ap- 
proach can serve as a good starting point for more 
sophisticated approaches correcting for errors in, e.g., 
the excitation spectrum, such as time-dependent DFT?^ 
the GW approximatiouf^i^i^S or self-interaction cor- 
rected DFT 4^1^ In weakly coupled molecular conduc- 
tors electron-electron interaction effects play a significant 
role. While some Coulomb blockade effects have been 
described using spin-density functional theory,— the cor- 
relation effects are more complicated to treat. In this di- 
rection the addition of a Hubbard-like term on top of the 
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DFT Hamiltonian has been used4^ These more advanced 
developments often come at the price of hmitations to 
the size of the systems that feasibly can be handled. It 
is therefore interesting to investigate to what extent the 
conventional DFT-NEGF can be used to model various 
transport properties. 

In this paper we present a scheme for including the 
effects of e-ph interaction into one such DFT-NEGF 
method for electronic transport. Specifically, we describe 
in detail our implementation of methods based on a com- 
bination of the SiestaM and the TransiestaS^ DFT 
computer codes. Siesta provides the fundamental imple- 
mentation of Kohn-Sham DFT in an atomic basis set for 
systems described in a supercell representation (periodic 
boundary conditions). Transiesta, on the other hand, 
uses the Siesta framework to solve self-consistently the 
Kohn-Sham DFT equations for the nonequilibrium elec- 
tron density in the presence of a current flow, taking 
into account the full atomistic structure of both device 
and electrodes (no periodicity in the transport direction). 
We describe how the Siesta and Transiesta methods 
have been extended for inelastic transport analysis, which 
involves the calculation of (i) relaxed geometries, (ii) vi- 
brational frequencies, (iii) e-ph couplings, and (iv) inelas- 
tic current-voltage characteristics up to the level of the 
self-consistent Born approximation (SCBA). We also de- 
scribe approximations leading to a lowest order expansion 
(LOE) of the SCBA expressions, which vastly simplifies 
the computational burdeni^i^ 

While there have already been many studies de- 
voted to transport with e-ph interaction based 
on model Hamiltonians emphasizing various aspects 
of the trancport ,'^^-'^''-^^-'^^-^"-^-^-^^-^'^'^^-^^'^^-^'^'^^-^^-^"-^-^ 
there has only been a handful based on a complete first- 
principles description of all aspects of the e-ph transport 
problem (described below) . By this distinction wc intend 
to emphasize approaches where structural, vibrational, 
and transport properties are derived from the knowl- 
edge of the elemental constituents only, i.e., without any 
system-dependent adjustable parameters. So far these 
have almost entirely been based on DFT for the elec- 
tronic structure. 

In the tunneling regime the atomic resolution of the 
STM has been used to investigate spatial variations of the 
inelastic tunneling process through adsorbed molecules 
on metallic surfaces. Corresponding inelastic STM im- 
ages were simulated theoretically by Lorcnte and Pers- 
son with DFT and the Tersoff-Hamann approachi^^iS 
Also controlled conformational changes, molecular mo- 
tion, and surface chemistry induced by the inelastic tun- 
nel current in STM have been addressed i^^j^i^ 

More recently the regime where an atomic-scale con- 
ductor is more strongly coupled to both electrodes 
has also been investigated. Based on a self-consistent 
tight-binding procedure with parameters obtained from 
DFT^ Pccchia et al. considered vibrational effects in 
octancthiols bonded to gold electrodes using NEGF and 
the Born approximation (BA) for the e-ph interaction!^ 



Solomon et al. further used this method to simulate the 
experimental lETS spectra of Wang et aL ^^^^^ Sergueev 
et al. studied a 1,4-benzenedithiolate molecule contacted 
by two aluminum leads This study addressed the bias 
dependence of the vibrational modes and e-ph couplings, 
but not the inelastic current itself. While the vibra- 
tional spectrum was found to be almost unchanged, a 
significant change in the e-ph couplings was found at 
high bias voltages (Vbias > 0.5 V). Chen et al. studied 
inelastic scattering and local heating in an atomic gold 
contact, a thiol-bonded benzene, and alkanethiolsi^SiZliiZ^ 
The inelastic signals were calculated using a golden-rule- 
type of expression and the DFT scattering states where 
calculated using jellium electrodes However, contrary 
to experiments and most calculations on molecules — for 
example Rcfs .20.,2I«67.468ii74ii75i.76i — they predict conduc- 
tance decreases by the phonons for alkanethiols. Jiang 
et al. used a related golden-rule approach for molecular 
systemsi^i Troisi et al. suggested a simplified approach 
from which lETS signals can be calculated approximately 
based on ab initio calculations for an isolated cluster and 
neglecting the electrodesi ^^i^^ This scheme was shown to 
be suitable for the off- resonance regime, i.e., when the 
molecular levels are far away from the Fermi level. Their 
results compare well with experiments by Kushmerick 
et ali^ During the development of the scheme presented 
here, we studied the same molecular systems with simi- 
lar resultsj^ii^ We also used it to model inelastic effects 
that can be observed in atomic gold wires 

The paper is organized as follows. In Sec.|TT]we commu- 
nicate our first-principles approach to obtain a Hamilto- 
nian description of a vibrating atomic-scale device bridg- 
ing two metallic contacts, such as schematically shown in 
Fig. [TJ Specifically we describe the use of Siesta to cal- 
culate vibrational modes and e-ph couplings. Section Hill 
addresses the NEGF formalism used to calculate the in- 
elastic electron transport in steady state as well as the 
SCBA and LOE schemes for the e-ph interaction. Elec- 
trode self-energies are obtained using the Transiesta 
scheme. Wc further discuss local heating effects and how 
various broadening mechanisms of the inelastic signal can 
be addressed. The main steps of the method presented 
in Sec. |TT] and IIIIl and how these depend on each other, 
are schematically clarified in Fig. [21 In Sec. IIVI and [Vl we 
illustrate our approach by corroborating and extending 
our previous studies of atomic gold wires and hydrocar- 
bon molecules. Section ITVl gives results for an extensive 
set of calculations for atomic gold wires of varying length 
and strain conditions. From these calculations wc iden- 
tify a number of physical effects, e.g., the evolution of a 
vibrational selection rule that becomes more pronounced 
the longer the wire is. Section |Vl illustrates that our 
method is applicable to a wide range of systems, here 
exemplified by different hydrocarbon molecules between 
gold surfaces. Both applications also underline the use- 
fulness of the LOE scheme, which wc validate by a com- 
parison the full SCBA calculation. Finally in Sec. IVII we 
provide a summary of the paper and an outlook. 
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FIG. 1: Schematic of two generic system setups, (a) To cal- 
culate vibrational frequencies and e-ph couplings with Siesta 
we use a supercell setup with periodic boundary conditions 
(BCs) in all directions. The cell contains the device region D 
and possibly some additional atom layers to come closer to a 
representation of bulk electrodes. The dynamic atoms are a 
relevant subset of the device atoms for which we determine 
the vibrations, (b) In the transport setup we apply the Tran- 
SIESTA scheme where the central region D is coupled to fully 
atomistic semi-infinite electrodes via self-energies, thereby re- 
moving periodicity along the transport direction (the periodic 
BCs are retained in the transverse plane). 



II. ELECTRONIC STRUCTURE METHODS 



In this section wc describe our first-principles method 
to obtain a Hamiltonian description of a vibrating 
atomic-scale device bridging to two metallic contacts. 
The framework is the Density Functional Theory (DFT) 
and its numerical implementation in the computer code 

SlESTA^i 



Vibrational Hamiltonian 



The physical situations which we typically want to de- 
scribe can schematically be represented as a central de- 
vice region D which is coupled to semi-infinite electrodes 
to the left (L) and right (i?) . This generic setup is shown 
in Fig. Mb). 

We assume that the whole system under consideration 



can be described by the following Hamiltonian 
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where cj and 6^ are the electron and phonon creation op- 
erators, respectively. Here is the single-particle mean- 
field Hamiltonian describing electrons moving in a static 
arrangement of the atomic nuclei, i/pj^ is the Hamilto- 
nian of free uncoupled phonons (oscillators), and -ffo-ph 
is the e-ph coupling within the harmonic approximation. 
For simplicity, we present in this paper a formulation 
for spin-independent problems. The generalization to in- 
clude spin-polarization is straightforward. 

The Hamiltonian Eq. ([1]) naturally arises from the 
adiabatic approximation of Born-Oppenheimer in which 
the timescales of electronic and vibrational dynamics 
are separated.— Since the electrons move on a much 
shorter timescale than the heavy nuclei, the adiabatic 
approximation states that the electronic Hamiltonian de- 
pends parametrically on the nuclear coordinates, i.e., 
that He = Hc{Q) where Q = R — R° is a displacement 
variable around the equilibrium configuration R*'. Next, 
limiting ourselves to small displacements we can expand 
the electronic Hamiltonian to lowest order in Q 
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where index / runs over all dynamic nuclei and v = x,y^z 
over spatial directions. Imposing a transformation into 
normal mode coordinates (and the usual canonical quan- 
tization of position and momentum operators) we can 
rewrite Eq. ^ into 
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where Mj is the mass of ion / and v''^ ~ {^lu} ionic 
displacement vector of normal mode A with frequency lox 
normalized according to v'^ • = 1. From Eq. ([3]) we 



identify the e-ph coupling matrix elements of Eq. (jldp as 



Mf, 



E( 

Iv 



2Miujx 



(4) 



In the following sections we describe how wc determine 
the detailed geometry, the vibrational modes, and the 
e-ph couplings from DFT. 
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FIG. 2: Flow diagram for the complete analysis of the inelastic transport properties of an atomic structure. 



B. SIESTA approach and geometry optimization 

In our numerical approach we use the Siesta imple- 
mentation of DFT4i This code treats exchange and cor- 
relation within the local density approximation (LDA) or 
the generalized gradient approximation (GGA) . The core 
electrons are described with pseudopotcntials. 

The main reason why Siesta is particularly suitable 
starting point for transport calculations is that the va- 
lence electrons are described in a localized basis set that 
allows for an unambiguous partitioning of the system into 
leads and device, cf. Fig. [Ijb), thereby making it possi- 
ble to calculate the flux of electrons (the necessity of 
this partitioning for transport calculations is discussed 
further in Sec. IIII|) . The basis orbitals {\i)} are strictly 
localized approximations to atomic orbitals with a given 
cutoff radius and centered at the positions of the nuclei of 
the structure. Importantly, this local electronic basis is 
nonorthogonal with overlap matrix elements Sij = 

In this tight-binding like basis we use the Kohn-Sham 
Hamiltonian from Siesta as the mean-field Hamiltonian 
in Eq. (jlbp . We initially construct a periodic superccU 
[Fig. [IJa)], and use it as an approximation to the full 
transport setup [Fig.[ljb)] for relaxing the device atoms, 
and to obtain vibrational frequencies and e-ph couplings. 
We note that this step leads to a determination of the 
quantities in equilibrium. In principle, these could also 
be calculated under nonequilibrium conditions by retain- 
ing the full transport structure of Fig. [Hb). Recently, 
Scrgucev et al. showed this to be important for relatively 
high voltages [eV ^ Hcux)^ However, for the low-bias 
regime considered in this paper the equilibrium calcula- 
tion is sufficient. 

A fairly accurate relaxation is an important prerequi- 
site for the subsequent calculation of vibrational modes. 



The atoms in the device region are therefore typically 
relaxed until the forces acting on the dynamic atoms all 
are smaller than Fi^{R°) < i^niax = 0.02 eV/A. Com- 
pared with other error sources in the calculations little is 
gained by lowering this criteria. 



C. Vibrational modes 

The starting point for our description of the nuclear 
vibrations is the Born-Oppenheimer total energy surface 
i?(R) (BOS) and its derivatives with respect to the nu- 
clear coordinates. For a thorough review on phonons 
from DFT we refer the reader to the paper by Baroni 
et alM- From the BOS we define the matrix of interatomic 
force constants (usually called the Hessian or dynamic 
matrix) as 



d^E{R) 



dRi^dR 



J^i 



(5) 



where R = {R/} denotes the full set of nuclear coordi- 
nates, and R/ = {Riu} the coordinates of nucleus / with 
mass Mj (not to be confused with the e-ph coupling ele- 
ments M^j). Within the harmonic approximation we can 
write the time-dependent displacement variable as 

Qj{t) = R/(i)-R? = Q/e'"*. (6) 

Inserting Eq. ([5]) and ^ into Newton's second law of 
motion 



we have 
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(8) 



5 



30 



> 

g 20 



u 

a 
o 

Oh 



-10 



Au dimer 

Pt dimer 



Acetylene Ethane 
Ethylene 



20 







SIESTA 
► Typical 
■ Accurate 
" ♦ Measured " 











400 
300 
200 
100 












ill 








\ ►♦ 








/ 


6' 


u 



FIG. 3: (Color online) Vibrational frequencies calculated for 
some simple molecules (Au2 and Pt2, acetylene C2H2, ethy- 
lene C2H4, and ethane C2H6). The results obtained directly 
from Siesta are shown together with those of our scheme 
(typical/accurate) based on the correction Eq. (|13[) . The dif- 
ferent calculational settings are described in the text. For 
comparison the experimentally measured values of the fre- 
quencies are also given.— To indicate the accuracy of 
the calculations the numerical values for the zero-frequency 
modes (translation/rotation) are included, where negative 
values correspond to imaginary frequencies. 



Introducing boldface notation also for matrices we can 
rewrite Eq. ([8| to the following ordinary eigenvalue prob- 
lem 



(w2i_W)v = 0, 



(9) 



where the mass-scaled matrix of interatomic force con- 
stants is 



(10) 



and v/ = v^M/Qj. Thus, the vibrational frequency lo\ 
and mode = {vj} belong to the eigensolution {uo\, v^) 
to Eq. ^ where we normalize the vectors as • = 1. 

Atomic forces F/ = {-F/i/} are directly obtained by 
Siesta along with the total energy calculation4^ This 
allows us to approximate the dynamic matrix by finite 
differences ( "frozen phonons" ) , either by 



7=r(±) 



Fi4±Qj^)-Fi,{0) 



or, numerically more accurately, by 

Fi,{Qj^)- Fi,{-Qj^) 



2Q 



(11) 



(12) 



where the bar denotes the finite difference approxima- 
tion. The quantities in Eq. (fTTI) and are thus readily 
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FIG. 4: (Color online) Convergence of calculated vibrational 
frequencies for a 4-atom Au wire with the most important 
DFT settings. For each of the two choices for the vibrational 
region (as indicated with boxes) the reference calculation — 
carried out with SZP, a 200 Ry real space grid energy cutoff, 
and 0.02 A finite displacements — and other three separate 
calculations (with one of the settings improved at a time) 
yield essentially the same results for the phonon energy hcux 
versus mode index A. 



determined. Typically we use a finite displacement of the 
dynamic atoms in each spatial direction oi Q = ±0.02 

A. 

While the Siesta calculations for Ciu-j,i are generally 
straightforward, we have observed that Siesta has dif- 
ficulties in estimating the change in force on the atom 
that is being displaced. This problem relates to the so- 
called egg-box effect, i.e., the movement of basis orbitals 
(which follows the nuclear positions) with respect to the 
real space integration grid.— As a result, phonons cannot 
be accurately obtained directly from Ci^^j^. To circum- 
vent this technicality we impose momentum conservation 
(in each direction v) via ^Fj^, = 0, which then de- 
termines the diagonal elements according to 



I = J 



(13) 



where the ii'-sum runs over all atoms in the supercell. Fi- 
nally, since d'^E/ dRi^dRj^ = d'^ E / dRj^dRiu we apply 
a numerical symmctrization of the force constants in the 
dynamic region. As a check we always verify that the fre- 
quencies calculated from the dynamic matrices with for- 
ward, backward, and combined displacements [Eq. (fTTI) 
and ()12|) ] are roughly the same, indicating that the har- 
monic approximation is not violated with the given dis- 
placement amplitude Qj^. 

The eigenvalues {ujj^ corresponding to the symmetric 
matrix W are real numbers. Some of these may however 
become negative leading to imaginary frequencies {ci^a}, 
indicating that the atomic configuration is in fact not 
describing a true energy minimum of the BOS. We shall 
denote such imaginary phonon frequencies by negative 
values in Fig. [3] and [51 
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A comparison between calculated and experimen- 
tally measured vibrational frequencies for some simple 
molecules is shown in Fig. [3l Specifically we include 
both the frequencies obtained directly with Siesta (from 
C'/^iJ/i) as well as those of our scheme based on the cor- 
rection Eq. ([in)) . In the calculations for the dimers the 
important settings correspond to cither a 200 Ry cut- 
off for the real space grid integration and a singlc-zeta 
plus polarization (SZP) basis set (Siesta/ typical), or a 
400 Ry cutoff and a double-zeta plus polarization (DZP) 
basis set (accurate). For the hydrocarbon molecules the 
settings are 200 Ry cutoff and DZP basis set. In all cal- 
culations the displacement amplitude is Qj^ = 0.02 A. 
The figure illustrates that our scheme presented above 
leads to a quite accurate description of the vibrational 
frequencies. We thus see no need to resort to a fre- 
quency scaling which is sometimes invoked in DFT cal- 
culations. Further, the figure shows that the use of 
momentum conservation for correcting elements in the 
Siesta dynamic matrix improves the calculation, in par- 
ticular the determination of low frequency modes (in- 
cluding the zero- frequency rotation/translation modes of 
isolated molecules). 

As an illustration of the convergence of the phonon 
energies with respect to some important DFT settings for 
larger systems, we show in Fig. 2] the calculated phonon 
energies for two different sizes of the dynamic region of 
a four atom gold wire (shown in the insets). We obtain 
almost identical frequencies by increasing the real space 
integration grid cutoff from 200 Ry to 300 Ry, by using a 
DZP basis set instead of a SZP, or by changing the finite 
displacements Qj^ from 0.02 A to 0.01 A. We expect the 
overall accuracy of these calculations to be representative 
not only for isolated molecules but also for larger periodic 
systems as well as systems involving other elements. 



-J2{r\k)iS-%i{l\H,\j) 

kl 

-Y,(^\H,\k){S~'),i{l\f). (16) 

kl 

The first term on the right hand side in Eq. can be 
approximated by finite differences of Hamiltonian matri- 
ces. The factors {i'\k) and are derivatives of the 
orbital overlaps, which we readily determine from finite 
differences via six separate runs that include both the 
original structure as well as the whole structure displaced 
by ±Qj^ along each spatial direction. We note that with 
the calculation of {i'\k) and we avoid the further ap- 
proximations for the e-ph couplings proposed by Head- 
Gordon and TuUy^ which we have used previously!^ 

In some cases, if one works with a relatively small su- 
percell, the calculated Fermi energy may change slightly 
between the displaced configurations of a given system. 
Since the real physical systems are essentially infinite, 
such shifts in the Fermi energy are artificial finite-size 
effects. To compensate for this we choose to measure 
all energies with respect to the Fermi energy of the re- 
laxed structure = eF(R-°), i-e., to shift the displaced 
Hamiltonians according to 



H(g/,) = H(g,,)-[eF(Q/.)-4]s(Q/.). 



(17) 



The finite difference approximation to the first term in 
Eq. (fT6|) — the derivative of the Hamiltonian matrix — 
may thus be written as 



D. Electron-phonon couplings 



In order to compute the e-ph coupling matrices M = 
{{M^j}} we have modified Siesta to output the Kohn- 

Sham Hamiltonian matrices H(Q) = {{(ijifob)}} for 
each of the displaced configurations. 

The complicated part of the e-ph couplings in Eq. (g]) 
is the evaluation of matrix elements of gradients of the 
Hamiltonian operator. Rewriting this part as 

(^i^b) = (14) 

where = d\i)/dQii, represents the change in basis 
orbitals with displacements, and using the identity 

Y,\^){^-%{J\^l, (15) 



where S = is the overlap matrix, we arrive at a 

form suitable for numerical evaluation 



9H 



« ^|H(Q,,)-H(-g,,) 

Q=0 ^WIu ^ 

-[eF(Q7.)-eF(-Q/,.)]S"}, 



(18) 



thereby completing the necessary steps to evaluate the 
e-ph coupling matrix elements. 



III. ELASTIC AND INELASTIC TRANSPORT: 
THE NEGF FORMALISM 

In this section we describe how the NEGF formal- 
ism is used to calculate the stationary electron trans- 
port through a region in space with an e-ph interaction. 
The basic ideas go back to the seminal work by Garoli 
et ali^ but we shall use the later formulation by Meir 
and Wingreen.^iMiSS 

The starting point in the NEGF approach is a for- 
mal partitioning of the system into a central device re- 
gion (where interactions may exist) and nonintcracting 
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leadsii2£ This partitioning was sketched in Fig.[Tfb). The 
e-ph interaction is treated with diagrammatic perturba- 
tion theory. Below we describe the SCBA as weh as 
further approximations leading to the computationally 
inexpensive LOE scheme. In addition, we discuss local 
heating effects and how various broadening mechanisms 
of the inelastic signal are addressed. 



A. System partitioning 

The physical system of interest sketched in Fig. [TJb) is 
infinite and non-periodic. For this setup let us initially 
consider the electronic and vibronic problems separately 
and return later to the treatment of their mutual inter- 
action. 

The use of a local basis in Siesta allows us to parti- 
tion the (bare) electronic Hamiltonian H = {{-ffj"}} and 
overlap matrix S = {{5,/,}} into 



Hl Hld 
H = I Hdl Hd H-dr 
Urd Hr 

Sl Sld 
Sl>l Sd ^dr 
SjiD Sji 



(19) 
(20) 



in which the direct couplings and overlaps between leads 
L and R are strictly zero (provided that the central region 
is sufficiently large). 

In a similar fashion, since interatomic forces are short 
ranged, the mass scaled dynamic matrix W [Eq. ^] can 
be partitioned into 



Wi Wld 
W = I Wdl Wdr 

Whd Wb. 



(21) 



where the direct coupling between leads L and R is ne- 
glected. 

The infinite dimensionality of the electronic and vibra- 
tional problem can effectively be addressed with the use 
of Green's function techniques. For the electronic part 
we define the retarded electronic single-particle Green's 
function G'''''(e) as the inverse of [{e + 177)8 — H] where 
77 = 0"*". It is then possible to write its representation in 
the device region D as 



[(£ + z,7)Sd Hd S2(e) - S^^X^r', (22) 



where the self-energy due to the coupling to the left lead 
is Si(e) = (Hdl - eSDL)gl{£){JiLD - £Sld) and sim- 
ilarly for the right lead. Here, g]^{e) is the retarded elec- 
tronic "surface" Green's function of lead a = L,R which 
can be calculated effectively for periodic structures by 
recursive techniques.— The quantities SJj(e) are directly 
available from TransiestAi^ Note that Green's func- 
tions calculated without the e-ph interaction are denoted 
with a superscript "0" . 



Similarly, for the vibrational part we can define the 
retarded phonon Green's function D°'''(tt') as the inverse 
of [{uj + irf)^! — W], and write its representation in the 
device region D as 



D"/(c.) = [(c. + z,?)2i-wc-n£(c.)-n5,i 



(23) 



where the self-energies due to the coupling to the left 
and right regions are — 'W oL^^Li'-^y^ ld and 

n^(w) = WD7^d^(w)WflX), respectively. Here, d^(u;) 
is the retarded phonon "surface" Green's function which 
again can be calculated by the recursion techniques men- 
tioned above. 

Note that the boldface matrix notation used for both 
electronic and vibrational quantities refers to different 
vector spaces: Indices in the electronic case refer to the 
basis orbitals and in the phonon case to real space co- 
ordinates. In addition, the electronic problem is treated 
directly in a nonorthogonal basis. The validity of the 
nonorthogonal formulation has be en disc ussed for the 
elastic scattering problem in Refs. 
cently including interactions in Ref. 

Since we are interested in the interaction of the 
electronic current with vibrations localized in the de- 
vice region, we invoke the ansatz that — to a first 
approximation — we can disregard the phonon lead self- 
energies IIJj(a;) and only describe the device region by 



87l88l and more re- 



D°/(u.) 



which in terms of the eigensolutions (lj^ , v 
can be written in a spectral representation 



(24) 

■^^ to Eq. (HI) 



where the free phonon Green's functions are^i 

1 1 



UJ ~ ujx ± if] ui + uj\ ± ir] 
+ i{nx) + l)S{uj±ujx)], 



(25) 

(26) 
(27) 



with (nx) being the expectation value of the occupation 
in mode A. The lesser and greater Green's functions 
stated above are used in Sec. IIII Dl (transformed into en- 
ergy domain via lo fuv). 

The validity of the approximation Eq. (|24[) can be 
investigated by calculating the correct phonon Green's 
function according to Eq. and then project the cor- 
responding local density of states (per energy via 1-^ e) 
onto each eigenmode v** of the dynamic region (with fixed 
electrodes), i.e., to determine 



Bxis) = -4e3m[(v^)^D237e)v^ 



(28) 
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satisfying the sum rule 



Jo 



(29) 



If the mode is a true locahzed modes for the extended 
system, then the projection B\{e) resembles a sharp res- 
onance around the phonon energy hiu\. In practice, {v^} 
are not exact eigenmodes of the extended system, and the 
resonances hence acquire finite widths. This broadening 
characterizes the damping (within the harmonic approx- 
imation) of the modes by the coupling to the electrodes. 
If the broadening is small compared with the phonon en- 
ergy (weak coupling to the bulk), then the projection can 
be described by a Lorentzian 



Bx{e) 



damp 



[e - fuvx 



amp 



(30) 



where fi-7damp ^''^^^ width half max (HWHM) value 

that transforms in time domain into an exponential decay 
of the phonon population with an average lifetime t^^^ = 

1 /7damp- ^ill return to the question of a finite phonon 
lifetime in Sec. [TnFl and HVEl 



B. Calculation of the current 

Our transport calculations are based on NEGF 
techniques and in particular the Meir-Wingreen 
formulation.— The steady-state (spin-degenerate) 
electrical current and the power transfer Pq, to the 
device from lead a ^ L, R can generally be expressed as 



2e(A^„) 



-2e 



Pa 



— oo 
oo 



de 
2^ 



2 f°° de 
Tr[S<(e)G>(e)-S>(e)G<(£)], 



(31) 

(32) 
(33) 



where Na is the electronic particle number operator of 
lead a, G^(e) the full lesser (greater) Green's function in 
the device region D (including all relevant interactions), 
and (e) the lesser (greater) self-energy that represents 
the rate of electrons scattering into (out of) the states 
in the device region D. We assume that the leads are 
unaffected by the nonequilibrium conditions in the device 
(this may be tested by increasing the device region). We 
can then use the fluctuation-dissipation theorem to write 
the lead self-energies as^i 



Sf(e) 



inF{e - Ha)Ta{e) 
i[nF{e ~ Ha) - l]r„(e) ' 



(34) 



where nF(e) = l/[exp(e/kBr) -I- 1] is the Fermi-Dirac 
distribution, /Xq, the chemical potential of lead a, and 

r„(e) ^ ^[^l{e) - S^(e)] = z[S>(e) - S<(e)], (35) 



describes the broadening of the device states by the cou- 
pling to the lead. 

The lesser and greater Green's functions are gener- 
ally related to the retarded and advanced ones via the 
Keldysh equation 



G>(£) = G2,(e)S>,(e)G^(e), 



(36) 



where S^j(e) is the sum of all self-energy contribu- 
tions (leads, interactions, etc.). Further, in steady- 
state situations time reversal symmetry relates the ad- 
vanced Green's function to the retarded one via Gfj(£) = 



C. Elastic transport 

If we consider a two-terminal setup with no interac- 
tions in the device region D, then the current expression 
simply reduces to the Landauer-Biittiker formula where 
Eq. ([33l) becomes 



taie) = [npis - ^J.L) - npie - hr)] 

X Tr[ri(e)G°/(e)rH(e)G2,'^(£)]. (37) 

Transiesta allows one to calculate the transmission 
function t{e) = ti^{e) — ti?(e) under finite bias condi- 
tions, i.e., with an electrostatic voltage drop over the 
device and different chemical potentials of the two leads. 
Due to the electrostatic self-consistency, this implies that 
the lead self-energies, e.g., H^{e), and Hamiltonian H 
depend parametrically on the external bias voltage V. 
These charging and polarization effects caused by the 
electrostatic voltage drop^i are fully treated in Tran- 
siesta at finite bias. Although it is relatively straight- 
forward to include these effects, it is computationally 
demanding for the inelastic calculation presented below. 
We have therefore neglected the voltage dependence and 
used the zero-bias self-energies and Hamiltonian in our 
inelastic calculations in the low-bias regime. In the case 
of metallic leads and a small applied bias (of the order 
of vibrational energies) we expect this approximation to 
be accurate. However, sufficiently large biases have been 
shown to influence the atomic structure^ as well as the 
e-ph couplings 



D. Self-consistent Born approximation 

Let us turn to the problem of the e-ph coupling. In 
order to use Eq. ([3T|) and ([32|) we need the full Green's 
functions Gj3(e) taking the e-ph interaction into account. 
Our approach is the SCBA where the phonon self-energy 
to the electronic system is described by the diagrams 
shown in Fig. [S}^ We note that in this work we ignore 
the phonon renormalization (pair bubble diagram) by the 
e-ph coupling. 
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FIG. 5: The lowest order diagrams for the phonon self- 
energies to the electronic description. The "Hartree" (a) 
and "Fock" (b) diagrams dress the electron Green's functions 
(double plain lines). The phonon Green's functions (single 
wiggly lines) are assumed to be described by the unperturbed 
ones, i.e., we ignore the e-ph renormalization of the phonon 
system. 

We write the phonon self-energies from mode A as^i^S 

(38) 

S;i,A(e) = ^[Sp\A(£)-Sp\A(e)] 

-^We'{Sp\A(£')-Sp\,(£')}(£), (39) 

where the retarded self-energy has been written in terms 
of the lesser and greater self-energies using the Kramers- 
Kronig relation He'{G'^(e')}(e) = iG^'ie). The func- 
tional 7i represents the Hilbert transform described in 
App.El 

The Hartree diagram Fig. [Slja) does not contribute to 
the lesser and greater phonon self-energies; this is because 
energy conservation implies that the wiggly line corre- 
sponds to a factor (i?(A,e' = 0) = Oi^ It does, however, 
lead to constant term for the retarded self-energy which 
can be understood as a static phonon-induced change in 
the mean-field electronic potentiah ^'^i^" From Eq. ([55]) 
we note that our retarded self-energy has the limiting 
behavior lime^oo Spj^ _;^(e) = 0. This is also the limits of 
the Fock diagram Fig. [5{b) if one calculates it directly 
with the Langreth rulesi ^^i^" We therefore conclude that 
Eq. ([M)) gives exactly the Fock diagram. Ignoring the 
Hartree term is reasonable since its small static potential 
shift might be screened (at least partially) if it had been 
included on the level of the DFT self-consistency loop. 
Further, the Hartree diagram docs not lead to a signal at 
the phonon threshold voltage. 

The full device Green's functions G^>(e) are related 
to G°/(e), SS^(£), and S;,^(e) ^ Ea ^txi^) via the 
Dyson and Keldysh equations^ 

GL(£) = G^'^(£) + G"/(£)S;;h(£)GL(£), (40) 
G%{e) = G^,(s)[sf(£) + s|(e) + s|(e)]G^(e). (41) 

The coupled nonlinear Eqs. p8 |) -(|4T |) have to be solved 
iteratively subject to some constraint on the mode popu- 
lation {nx) appearing in dg'(A,e), cf. Eq. (^7)) . For weak 
e-ph coupling we thus approximate the mode occupation 
{n\) by the steady-state solution to a rate equation de- 



scribing the heating of the device 

("a) = ^y- -7damp[('^A) -nB(fttJA)], (42) 

where nB(e) = l/[exp(£/kBr) — 1] is the Bose-Einstein 
distribution, px the power dissipated into mode A by the 
electrons, and 7^amp ~ ^/'''ph ^ damping parameter re- 
lated to the average lifetime of the phonon, e.g., by cou- 
pling to bulk vibrations. 

In steady state the power transferred by electrons from 
the leads into to the device must balance the power trans- 
ferred from the device electrons to the phonons, i.e., 

A 

From the particle conservation condition^S 

Tr[S<,(e)G>(e) - S>,(£)G<(£)] = 0, (44) 
we can define the quantity px as 

xTr[S<h,,(£)G>(£) - S>, ,(e)G<(e)], 

which consequently obeys Eq. (|43)) . Wc note that in this 
way we basically define 3A^ quantities from a single equa- 
tion for X^Af^A only; different definitions could in prin- 
ciple also fulfill the power balance. However, to lowest 
order in the e-ph coupling our definition Eq. (|45p is un- 
ambiguously the power transferred to mode A. 

From Eq. (|^^ we can identify two regimes, (i) the ex- 
ternally damped limit (7damp much larger than electron- 
hole (e-h) pair damping 7c_h) where the populations 
arc fixed according to the Bose-Einstein distribution 
{nx) = UBifiuJx)-, and (ii) the externally undamped limit 
(7damp = and hence from Eq. (|^ that px = 0) where 
the populations vary with bias such that no power is dis- 
sipated in the device, i.e., Pl + Pr = 0. It is instructive 
to note that px includes both phonon emission and ab- 
sorption processes, which is the reason why a steady-state 
solution always exists. 

A typical situation that come close to the externally 
undamped limit is when the device vibrations fall out- 
side the phonon band of the bulk electrodes, i.e., when 
there is a significant mass difference between the device 
atoms and the electrode atoms. In this case the vibra- 
tions cannot couple directly (resonantly) to the bulk, and 
the damping, e.g., by anharmonic means, is likely to be 
much smaller than the coupling to the electrons. One 
important example is the hydrogen molecule clamped be- 
tween platinum contactsJ^ii^ 

To solve the SCBA equations Eqs. (l38l)-(l42l), we have 
developed an implementation in the programming lan- 
guage Python where the Green's functions and self- 
energies are sampled on a finite energy grid. The main 
technical challenges are discussed in App. |B] Finally we 
note that with the phonon self-energies Eqs. (|38p - ((55|) 
the current is conserved. This can be proven using the 
identity Eq. 
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E. Lowest order expansion 

The solution of the SCBA equations is a daunting nu- 
merical task for systems consisting of more than a hand- 
ful of atoms. However, for systems where the e-ph cou- 
pling is weak and the density of states (DOS) varies 
slowly with energy, we have previously derived the LOE 
approximation.— Here we elaborate on these results. 

The main computational burden of the SCBA origi- 
nates from the numerical integration over energy needed 
in the evaluation of the current and power expressions 
Eqs. ([3T|) -([32 |) . The LOE approximation assumes that 
the retarded and advanced single-particle Green's func- 
tions G^'^^" and lead self-energies TfJ"" are energy inde- 
pendent. We can then expand the current and power 
expressions to the lowest order (second) in e-ph cou- 
plings M"^ and perform the energy integrations analyti- 
cally. These integrals consist of products of Fermi-Dirac 
functions and their Hilbert transforms. The LOE thus 
retains the Pauli exclusion principle for fermionic parti- 
cles, which is necessary to model the blocking of phonon 
emission processes at low bias. 

In the LOE approximation, the total power dissipated 
into the phonon system P^OE ^ _|_ ^^^^ after 
lengthy derivations, be written as^ 



3LOE 



LOE 



(46) 



pY''' = hi^x {[nB(fi^A) - {nx)hth + -/LiV.T)} , (47) 



7c-h 



7^ 

I cm 



— ^Tr [M^AM^Al 



huJx [cosh( j^^) — 1] coth( 



2kBT) 



(48) 

.ysinh(^) 



^ft[cosh(|^)-cosh(^)] 

xTr [M^AlM^At^J , (49) 

where the Bose-Einstein distribution nB(e) appears in 
Eq. (|T7|) due to the integration of Fermi-Dirac func- 
tions describing the electrons in the contacts. Here 
G = G^''(£f), r„ = r„(eF), and A = i{G - G^) are 
the noninteracting retarded Green's function, the broad- 
ening by contact a = L, R, and the spectral function at 
Ef, respectively. For convenience we have also defined 
the quantities A^ = Gr^Gt such that A = A^ -I- A^. 

The first term in Eq. (j47p describes the equilibrium 
energy exchange between the vibrational and electronic 
degrees of freedom (e-h pair damping 7g_jj of the vibra- 
tions); it tend to drive the phonon system towards the 
Bose-Einstein distribution. The second term appears in 
nonequilibrium and is related to an effective emission rate 
7gjjj of vibrational quanta under finite bias. At low tem- 
peratures (keT 0) this rate is given as 



7cm 



\eV\ - hwx 



\eV\ - hiux)TT [M^AlM^A^] , 

(50) 



where 0{x) is the step function; i.e., the net emission of 
phonons above the threshold grows linearly with the bias 




...Hi 



Left 



Right 



1 



d) 



FIG. 6: Schematic representation of the energy phase space 
available for scattering processes due to the Pauli principle. 
Phonon emission (a) and absorption (b) between scattering 
states originating from the left and right contacts. Figs, (c) 
and (d) correspond to phonon absorption between scattering 
states in the same contact. 



voltage. Furthermore, since Tr [M^AqM'^'A^] > 0, we 
find that 



Tr [M^AM^A] > 2Tr [M^AlM^A^] 



(51) 



We can use this inequality to derive an upper bound on 
the phonon occupation by solving the steady-state con- 
dition pY''^ = [cf. Eq. (|42|) with no external damping] . 
It simply bccomes2ii2^ 



(52) 



, , 1 \eV\ — fuox ^ , 

("a) < ^ ' 6{\eV\~hujx). 



To provide an intuitive understanding of Eqs. 
()52p consider the following arguments: The energy phase 
space available for phonon emission and absorption pro- 
cesses is limited by the Pauli principle, as sketched in 
Fig. [6l We divide the electronic phase space in two, cor- 
responding to scattering states incoming from either the 
left or the right contact. Without e-ph scattering these 
states are assumed to be populated up to the Fermi level 
£f (we take iil > IJ'R + and keT 0). Within 
this picture phonon emission can only take place from 
a populated state originating in the left contact to an 
empty state originating in the right contact, see Fig.lHlJa). 
Similarly, phonon absorption can be described by three 
different processes sketched in Fig. [n^b)-(d), again cor- 
responding to scattering from populated initial states to 
empty final states. 

The scattering rates for these processes are propor- 
tional to the energy window in which they can take place. 
Denoting the scattering rate per energy as djaa' /de, 
where a ~ L, R {a' = L, R) indicates the propagation 
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direction of the initial (final) scattering state, we can 
write the spontaneous plus stimulated emission power as 
fi^\{{nx) + 1) {eV — huj\) djLji/de and the ab- 
hujx{nx)[{eV + hhJx) djLR/de 
+liujx (d'jLL/de + djuu/de)]. The net power transfer 
from the electronic system to the phonon mode A is there- 
fore 



„LOE 

sorption power as pYi^ 



LOE LOE 

Pa, cm ^ P\,ah 



-{hwxY{nx) 



LL 



+hujx {eV — hujx) 



de de 



djRR 

de 



(53) 



A comparison with Eq. (|47p reveals that the term propor- 



tional to the occupation (nx) is bias independent due to a 
cancelation of phonon absorption by stimulated emission. 
Furthermore, the upper bound in Eq. ([52|) is directly mo- 
tivated by equating Eq. (|53p to zero (steady state) and by 
ignoring scattering processes with initial and final states 
propagating in the same direction (d'^aa/de). In addi- 
tion, a steady-state solution to Eq. (|42|) always exists 
because the phonon emission rate is always smaller than 
the total phonon absorption rate, and that emission pro- 
cesses are restricted to a smaller energy window than 
absorption processes. 



The LOE approximation, which above was applied to 
the power, also allows us to write the current through 
the device 1^°^ 



rLOE 



T-sym 
-^A 



GoVTr [Gr^Gtri] 



Gtr^G \ M^A^M^ + - (rj^G^M-^AM-^ - h.c. 



Y,rr{V.T,{nx))T, 

A 

^Xf '""(F, T) Tr [G^r^G {r^G^M^ [Ar - A^) + h.c.}] 

A 

hjjx + eV \ 



1 e ''B^ 



OO 

6 f diE 

^ / -^[nF{.!^)-'nF{e- eV)]He'{nY{e' + hwx)-nY{,e' -hujx)}{£), 



(54) 
(55) 
(56) 



r 



where the bias is defined via eV ~ ^r — ^il, and Go = 
2e?lh is the spin-degenerate conductance quantum. This 
expression is current conserving, i.e., calculating the cur- 
rent at the left and right contacts give the same result. 

The LOE expression for the current Eq. ([51)1 con- 
tains three terms, (i) the Landauer-Biittikcr term cor- 
responding to the elastic conductance, (ii) the "symmet- 
ric" term corresponding to symmetric conductance steps 
at the vibrational energies, and (iii) the "asymmetric" 
term corresponding to peaks and dips in the conduc- 
tance which are asymmetric with voltage inversion, see 
Fig. [71 For geometrically symmetric junctions, it can be 
shown that the asymmetric term vanishes exactly. Even 
for geometrically asymmetric systems we typically find 
that it is a very small contribution compared with the 
symmetric term. Furthermore, the sign of the conduc- 
tance step for the symmetric term in general shows an 
increase (decrease) in the conductance for low (high) con- 
ducting systems, e.g., vibrations usually help electrons 
through molecules while they backscatter electrons in 
atomic wires. This is discussed further for a one-level 
model in Ref. [ol. 

The LOE approximation is computationally simple 



and can be applied to systems of considerable size. Al- 
though the approximation is not strictly valid for systems 
with energy-dependent DOS, comparison with the full 
SCBA calculations shows good agreement even for sys- 
tems that have a slowly varying DOS (on the scale of vi- 
brational energies), e.g., the organic molecules connected 
to gold electrodes described below in Sec. |Vl The LOE 
approximation will certainly fail when sharp resonances 
(compared to the vibrational energies) are present within 
the order of phonon energies of the Fermi energy. How- 
ever, in this case Coulomb blockade physics is also ex- 
pected, which thus makes any DFT mean-field approach 
(including ours) questionable. 



F. Broadening mechanisms 

The width of the experimentally measured phonon sig- 
nal in the conductance is a combination of (at least) three 
broadening mechanisms, namely the intrinsic ones from a 
finite temperature and a finite phonon lifetime, as well as 
the one related to the modulation voltage used in lock- 
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FIG. 7: (Color online) Universal functions Eq. ([55|l and ([56|l 
giving symmetric and asymmetric phonon contributions to 
the conductance in the LOE, respectively. The differential 
conductance dl/dV and the second derivative d?I/dV^ are 
shown (in arbitrary units) for one phonon mode for three 
different temperatures (a) \i^T/hujx = 0.02, (b) \i-BT/Tujx = 
0.06, and (c) ksT/fic^A = 0.10. 



in measurements (to improve the signal-to-noise ratio) 
of the second derivative of the current with respect to 
the bias. These contributions do not add up triviaUy. 
However, as we show below, one can provide estimates 
for each of the different contributions which thus help to 
understand what effect is the dominant one. 

As can be seen in Fig. [71 the electronic temperature 
gives rise to a broadening of the vibrational signal. From 
Eq. dSSl) the full width half max (FWHM) in the second 
derivative of the current can be shown to be approxi- 
mately 5.4 k^T^M^ 

The effects of a finite phonon lifetime t\ — 1/ 7damp 



to a first approximation described by a convolution of the 
free phonon Green's functions with a Lorcntzian with a 
HWHM width of fi.7damp- Consequently, this convolution 
propagates to the phonon self-energies Eq. ((38)) and to the 
inelastic LOE corrections to the current, cf. Eq. and 
([55)1 . The FWHM broadening in the second derivative 
of the current is thus 2?i7damp- The intrinsic lincwidth 
of the phonon signal has also been discussed in a simple 
SCBA model by Galperin et al^ 

The broadening from the lock-in technique for mea- 
surements of the first or second derivatives of the current 
can be estimated in the following way. With a small har- 
monic modulation signal (with amplitude A = K-ms) 
applied on top of the bias voltage one can measure deriva- 
tives of the current. As shown in App. [C] the FWHM 
width induced by the lock-in measurement technique is 
2.45 Vrms and 1.72 Vi-ms for the first and second derivatives 




FIG. 8: (Color online) Generic gold wire supercells containing 
3 to 7 atoms bridging pyramidal bases connected to stacked 
Au(lOO) layers. As indicated on the figure, the electrode sep- 
aration L is defined as the distance between the plane in each 
electrode containing the second-outermost Au(lOO) layer. 



of the current, respectively (neglecting intrinsic broaden- 
ing). In other words, if d?I/dV^ is a (5-function, the 
experimentally measured FWHM width will be either 
2.45 Vrms or 1.72 Vmis, depending on whether the lock- 
in measurement is on the first or second harmonic. 



IV. ATOMIC GOLD WIRES 

Since the discovery in the late 1990s that gold can 
form free-standing wires of single atom a^°°'''^°"'^i^°^'''^°^ 
the mechanical, chemical, and electrical proper- 
ties of these atomic-scale systems have been extensively 

„-|-^j^^^ 15.48.78.104.105.106.107.108.109.110.111.112.113.114.115.116.117.118.119 .1 

For this reason we illustrate in this section our method 
described in Sec. |lT] and |TTT] by applying it to model 
inelastic scattering in atomic gold wires. We compare 
directly the results of our theoretical developments with 
the high-quality experimental data by Agrai't and co- 
workcrs .^^'^^^ They used a cryogenic STM to first create 
an atomic gold wire between the tip and the substrate 
surface, and then to measure the conductance against 
the displacement of the tip. From the length of the 
observed conductance plateau around Go — the signature 
that an atomic wire has been formed — it was possible 
to determine the approximate size as well as the level of 
strain of the created wire. Under these conditions Agrai't 
et al. then used point-contact spectroscopy to show that 
the conductance of an atomic gold wire decreases a 
few percent around a particular tip-substrate voltage 
(symmetric around zero bias) presumably coinciding 
with the natural frequency of a certain vibrational mode 
of the wire. With this inelastic spectroscopy method 
they could further characterize the conductance drop as 
a function of wire length and strain. 

To simulate these experiments, we study wires contain- 
ing different number of atoms and under varying stretch- 
ing conditions. The generic supercells used in the Siesta 
calculations are illustrated in Fig. [8] and consist of 3 
to 7 gold atoms bridging pyramidal bases connected to 
stacked Au(lOO) layers. We use a 4 x 4 supcrccU size in 
the plane transverse to the transport direction and define 
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the electrode separation L, as indicated on Fig. [H as the 
distance between the plane in each electrode containing 
the second-outermost Au(lOO) layer. The face-centered 
cubic (FCC) lattice constant for the bulk gold atoms is 
taken to be a = 4.18 AJ^ 

We generally use (unless otherwise specified) the 
Perdew-Burke-Ernzcrhof version of the GGA for the 
exchange-correlation functionalfi^ a split-valcncc single- 
C plus polarization (SZP) basis set with a confining 
energy of 0.01 Ry (nine orbitals corresponding to the 
5d and 6(s,p) states of the free Au atom), a cutoff 
energy of 200 Ry for the real space grid integration, 
and the F-point approximation for the sampling of the 
three-dimensional Brillouin zone. The interaction be- 
tween the valence electrons and the ionic cores are de- 
scribed by a standard norm-conserving TrouUier-Martins 
pseudopotentiaU^ generated from a relativistic atomic 
calculation (including core correction). We have found 
that these settings yield a reasonable compromise be- 
tween accuracy and computational cost. 



A. Geometry relaxation 

For a given electrode separation L the first calcula- 
tional step is to relax the geometry to obtain a local 
energy minimum configuration R'^. With the settings 
described above we relax both the outermost electrode 
layers, the pyramidal bases, and the wire atoms until 
all forces acting each of these atoms arc smaller than 
i^„,ax = 0.02 eV/A. 

Figure [9l^a) shows the relative differences in the Kohn- 
Sham total energy (cohesive energy) as the wires are elon- 
gated. We also show the numerical derivatives of these 
binding energy curves as a measure of the forces acting 
on the wire. The breaking force, defined as the energy 
slope of the last segment before breaking, is found be of 
the order 1 eV/A ^ 1.6 nN. This agrees well with the 
experimental results which have shown the break force 
for atomic gold wires to be close to 1.5 nN^iiiiiii^ 

In Fig. [5]Jb) we summarize the geometrical findings 
of the relaxation procedure by plotting the wire bond 
lengths and bond angles as a function of electrode separa- 
tion L. The figure shows that the short wires containing 
3 or 4 atoms adopt a linear structure over a wide range 
of electrode separations. The longer wires, on the other 
hand, are generally found to have a zigzag geometry only 
approaching a linear form when they are stretched close 
to the breaking point 

From the plot of the bond lengths between nearest 
neighbors in the wire one notices that the 4 and 6 atom 
wires have a more pronounced tendency to dimerizc than 
the wires with an odd number (due to left /right symme- 
try of the structures only wires with an even number of 
atoms should be able to dimerize). In three test calcula- 
tions with a 3 X 3 X 3 k-point sampling of the three dimen- 
sional Brillouin zone we generally achieve very similar 
atomic arrangements as compared to the F-point only. 
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FIG. 9: (Color online) Energetic, geometric, and conductive 
properties of atomic gold wires: (a) Kohn-Sham total energy 
(cohesive energy) vs. electrode separation, (b) bond angles 
and bond lengths, (c) phonon energies, and (d) elastic trans- 
mission at the Fermi energy calculated both for the F-point 
(colored open symbols) as well as with a 5 x 5 k-point sam- 
pling of the two-dimensional Brillouin zone perpendicular to 
the transport direction (black stars). 
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FIG. 10: (Color online) Generic transport setup in which a 
relaxed wire geometry — here a 7-atom wire with L — 29.20 
A — is coupled to semi-infinite electrodes. As indicated on the 
figure the vibrational region is taken to include the atoms in 
the pyramidal bases and the wire itself, whereas the device 
region (describing the e-ph couplings) includes also the out- 
ermost surface layers. 

However, these calculations, which are indicated with 
black crosses in Fig. Mi>), seem to reduce the dimeriza- 
tion tendency somewhat. 



B. Vibrational analysis 

We calculate the vibrational frequencies and modes as 
described in Sec. Ill CI With N vibrating atoms we thus 
find SN modes for a given structure. The phonon spec- 
trum for the wire is plotted in Fig. [9jc) , where negative 
values indicate modes with imaginary frequency implying 
the breaking of an unstable wire. The general trend is 
that the phonon energies diminish as the wires are elon- 
gated. This can be understood by considering that the 
effective "springs" between ions in the wires are softened 
as the bonds are stretched, which in turn result in lower 
energies. 

In the results to follow we generally take the wire and 
pyramidal base atoms as the dynamic region (as indi- 
cated in Fig. [To)) , i.e. these atoms are allowed to vibrate. 
For the 3- to 7-atom wires this leave us with 33 to 45 vi- 
brational modes. The corresponding e-ph couplings are 
calculated in a slightly larger device region containing 
also the outermost surface layer. This inclusion of an 
extra layer is necessary to represent the vibrational mod- 
ulation of the hopping between the pyramidal base atoms 
and the first surface layers. 



C. Elastic transmission 

In order to determine the transport properties of the 
wire geometries described above, we construct from the 
supercells shown in Fig. [5] new wire geometries which 
are coupled to semi-infinite electrodes as schematically 
illustrated in Fig. [Ijb). The resulting setup is shown in 
Fig. [To] for the case of a 7-atom long gold wire. As indi- 
cated on this figure we consider the device subspace to 
include the top-most surface layer, the pyramidal bases, 
and the wire itself. 



The elastic transmission evaluated at the Fermi energy 
Ef is calculated using Transiesta described in Ref. [2y. 
The results are shown in Fig. [H^d) both for the F-point 
(open symbols) as well as with a 5 x 5 k-point sampling 
of the two-dimensional Brillouin zone perpendicular to 
the transport direction (black star s). In corresp ondence 
with previous work, e.g., Refs. [87l ll04lll07lflTi we find 
that the total transmission is close to unity, except for 
the very stretched configurations where the transmission 
goes down somewhat. From Fig.[SJd) one observes a rea- 
sonable agreement between the F-point and the k-point 
sampled transmissions, particularly when the transmis- 
sion is close to one. Worst are the discrepancies for the 4 
and 6 atom wires, which also arc the cases where the 
transmission deviates most from unity. We subscribe 
these signatures to the so-called odd-even behavior in 
the conductance of metallic atomic wires, in which per- 
fect transmission is expected only for an odd number of 
atoms in a chain. For an even number of atoms the con- 
ductance should be lower Further, the observed dimer- 
ization is also expected to reduce the conductance (the 
Pcicrls instability for infinite metallic wires results in the 
opening of a band gap at the Fermi energy) . We also note 
that on an energy scale of the typical phonon energies the 
transmission function is to a very good approximation a 
constant around the Fermi energy. 



D. Inelastic transport 

Having determined the vibrational frequencies, the e- 
ph couplings, and the elastic transmission properties, we 
are in position to calculate the inelastic current as de- 
scribed in Sec. IIIIBI 

We start out by showing that the LOE and SCBA ap- 
proaches essentially predict the same inelastic signals for 
atomic gold wires, thereby reducing the computational 
expense in the detailed analysis to follow. For this pur- 
pose only we consider a computationally reduced problem 
where the device and dynamic atoms regions are mini- 
mized as compared with those generally adopted in this 
section. We will thus simply allow the wire atoms to vi- 
brate and take the device space as the wire plus pyrami- 
dal bases only. Compared with the electronic structure 
and phonon energies the thermal energy typically sets 
the smallest energy scale for variations in the Green's 
functions etc. Instead of using the experimentally rele- 
vant temperature of T = 4.2 K (or even less) we further 
simplify the calculations by taking T = 10.0 K for the 
moment since this requires fewer points on the energy 
grid, cf. App. [H 

The differential conductances as resulting from eval- 
uating Eq. ([^T|) with and without SCBA phonon self- 
energies as well as evaluating the LOE expression 
Eq. ((54)) are shown in Fig. [TT] The dotted curve is the 
purely clastic result (no phonon self-energy) and the cir- 
cles the full SCBA (including all vibrational modes in 
the externally damped limit 7damp ^ 7e-h of Sec. IIII D]) . 
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FIG. 11: (Color online) Elastic and inelastic differential con- 
ductance calculated at T = 10.0 K in a reduced device region 
for the 7-atom wire shown in Fig. 1101 The small variation in 
elastic conductance with bias (dotted curve) relates to a weak 
energy dependence of the elastic transmission function at the 
F-point around ep. The full SCBA calculation (circles) follows 
this trend and shows on top of it symmetric drops character- 
istic for phonon scattering. The LOE calculation (line) does 
not include the elastic variation but gives basically the same 
predictions for the inelastic signals as the SCBA with the elas- 
tic background signal subtracted (dashed curve). This illus- 
trates the agreement between the LOE and SCBA approaches 
for the inelastic contribution. 



The red line corresponds to the LOE. The elastic conduc- 
tance displays a slight variation with bias that relates 
to the weak energy dependence in the zero-bias trans- 
mission function at the F-point. The full SCBA calcu- 
lation clearly shows two symmetric conductance drops 
which are due to inelastic scattering against vibrations 
(we will return later to a discussion of the physics) . The 
LOE calculation does not include the elastic variation 
but gives basically the same predictions for the inelas- 
tic signals. This is clear from a comparison with the 
SCBA where the elastic background signal has been sub- 
tracted (dashed curve). Based on a number of such tests, 
and the fact that the e-ph couplings are weak (or more 
precisely, that the inelastic signal is a small change in 
conductance of the order 1-2 %), we conclude that the 
approximations leading to the LOE expressions are valid 
in the case of atomic gold wires. To appreciate this fact, 
we note that the SCBA curves in Fig. [TT] required ap- 
proximately 40 CPU-hours in a parallel job running on 4 
processors whereas the LOE results only required a few 
seconds on one processor. The LOE approach is thus 
justified for a full analysis of the 3- to 7-atom gold wires. 

Figure[T2]shows the calculated differential conductance 
of the 3- to 7-atom wires under different electrode sepa- 
rations L and in the externally damped limit. The de- 
vice region and dynamic atoms are here as indicated in 
Fig. [TOl and the temperature of the leads is T = 4.2 K. 
The curves display symmetric drops at voltages corre- 
sponding to particular phonon energies. The dominant 
inelastic signal moves towards lower energies and increase 



in magnitude as the wires are elongated. Furthermore, 
sometimes also a secondary feature is found below 5 meV, 
e.g.. Figs. [TT] and [T^ These observations are also charac- 
teristic for the experiments f^ii^ and in agreement with 
previous calculations 4^*2^ 

To extract the general trends on how the inelastic sig- 
nal depends on details in the atomic arrangement we 
present in Fig. [T3] our calculated data in different forms. 
In these plots we represent each phonon mode by a dot 
with an area proportional to the corresponding conduc- 
tance drop. The abscissa corresponds to the electrode 
separation whereas the ordinate is used to highlight cer- 
tain properties of the vibrational modes. In this way, 
Fig. fTHla) illustrates the mode frequency change with 
electrode separation. From a linear fit to the strongest 
signals we predict a frequency shift of —8.45 meV/A for 
the 5-atom wire falling off to —6.34 meV/A for the 7- 
atom long wire. Further, to understand the nature of 
the modes that infiucncc the electronic transport we can 
try to quantify some important characteristics. As it has 
previously been shown, longitudinal modes with an al- 
ternating bond length (ABL) character are expected to 
be the dominating onesFi^i^^ii^^ To measure the longi- 
tudinal part of a given vibrational mode we define 
a sum over z-components ^^/(^/z)^ — ^ where / runs 
over all dynamic atoms (the upper bound is due to the 
eigenmodes normalization v'^ • v'^ = 1). This quantity is 
shown in Fig. [TSFb) . The plot clearly expresses that the 
modes with the largest signals (large dot area) also have 
a strong longitudinal component. Further, to show that 
these modes also have ABL character, we also define a 
sum X]/>7 K/z ~ ^jzl where / and J are nearest neigh- 
bor atoms in the chain. This second quantity is shown 
in Fig. [T^T c). from which we learn that the important 
modes also have the largest ABL measure (the absolute 
scale is irrelevant). 

Another important aspect is whether the modes are 
really localized in the wire or not. Remember that our 
approach assumes that atoms outside the dynamic region 
are fixed. Therefore, if we have eigenvectors with signifi- 
cant amplitude near the boundary of the dynamic region, 
this assumption does not seem to be valid (most likely the 
eigenvector is not a true eigenvector of the real system). 
In other words, we want to make sure that the modes 
which arc responsible for the inelastic scattering are suf- 
ficiently localized "deep" inside the dynamic region. To 
show this we calculate J^i Vj • Vj < 1 where / runs over 
the 3 to 7 wire atoms. This quantity is represented in 
Fig. 113r d) and confirms that indeed the important modes 
are localized in the chain; particularly for the 5-, 6-, and 
7-atom wires the localization is almost perfect. 

In conclusion, from the results presented in Fig. 1131 we 
learn that the inelastic signal in the conductance is effec- 
tively described by a simple selection rule in which lon- 
gitudinal vibrational modes with ABL mode character — 
localized in the wire — are the main cause of the inelastic 
scattering. We are further able to quantify the frequency 
down-shift and signal increase with strain. 
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FIG. 12: (Color online) The differential conductance G and its derivative dG/dV calculated with the LOE approach for the 3- 
to 7-atom gold wires in the externally damped limit. The electrode separation L is indicated next to the conductance curves. 
As shown in Fig. [10] the device region includes the outermost electrode layer whereas the dynamic atoms are pyramidal bases 
plus wire. The temperature of the leads is T = 4.2 K. 



E. Vibrational lifetimes and local heating 

From Fig.fT^Td) wc get a hint about the damping of the 
modes from the couphng to bulk phonons. If a mode is 
localized "deep" inside the dynamic region this coupling 
is negligible and the mode is expected to have a long 
life-time, i.e., to be weakly damped by the coupling to 
the bulk. As discussed in Sec. IIII we can estimate this 
damping from the width of the phonon density of states 
projected onto the mode vector. 

As an illustration of this approach, we calculate the 
damping of the dominating ABL mode according to 
Eq. (|30p in the case of the 7-atom wire with electrode sep- 
aration L = 29.20 A. This mode, shown in Fig.fHTa). has 
a localization quantity (as defined above) of value 0.987, 
i.e., it is 98.7% localized in the wire. We begin by deter- 
mining the dynamic matrix of the whole wire supercell 
[Fig. [HKe)] as described in Sec. Ill CI To describe the bulk 
properties of gold we pick the intra-layer and inter-layer 
elements (inside the slab) in the dynamic matrix along 
the transport direction, and use recursive techniques to 
calculate bulk and surface phonon Green's functions. Be- 
cause of periodicity in the transverse plane — which gives 
rise to artificial sharp resonances in the spectrum — we 
broaden the phonon Green's functions by taking 77 = 1.0 
meV. This approach leads to the total phonon density 
of states (full black line) shown in the inset of Fig. [141 
This shape compares reasonably well with other calcula- 
tions and experimentSF^^ii^ The inset also shows the 
phonon density of states decomposed in the direction 



of the transport (dashed red curve) as well as in the 
transverse directions (dotted blue curve); the observed 
isotropy that is expected for bulk is actually quite satis- 
factory. Finally, we calculate the projected phonon den- 
sity of states Bx{lu) for the ABL mode of interest accord- 
ing to Eq. ([50]) . This projection on a discrete energy grid 
is shown in Fig.[T4|(open circles). By fitting a Lorentzian 
to the calculated data points we obtain a FWHM of 8 ficV 
and a shift in frequency by —6 /xeV. Based on these cal- 
culations we thus estimate the phonon damping to be of 
the order fi7damp = ^ fj,eV. In fact, this is rather a lower 
bound, since we have not included anharmonic contribu- 
tions etC;Si However, compared with the phonon energy 
we see that indeed 7damp "^a, and thus that the use of 
free phonon Green's functions in the SCBA self-energy 
Eq. ([SHU is justified. 

Let us next investigate the implications of a finite 
phonon lifetime on the local heating. This is done by 
solving the rate equation Eq. (|42p for the mode occupa- 
tion at a fixed bias voltage. For instance, the inelastic 
conductance characteristics (including heating) for our 
7-atom wire arc shown in Fig. [T5| for different values of 
the phonon damping 7damp (smooth colored lines). As 
seen in the figure, and as we have shown previously)^ 
the effect of the heating is to introduce a slope in the 
conductance beyond the phonon threshold voltage. This 
is because the nonequilibrium mode occupation increases 
the number of scattering events of the traversing elec- 
trons. Consequently the conductance goes down as the 
bias (and hence the occupation level) increases. The 
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FIG. 13: (Color online) Inelastic signals plotted as a function of the electrode separation L. Each mode is represented by a dot 
with an area proportional to the corresponding conductance drop. On the j/-axis we show (a) the phonon mode energy, (b) a 
measure of the longitudinal component of the mode, (c) a measure of the ABL character, and (d) a measure of the localization 
to the wire atoms only. The straight lines in plot (a) are linear interpolations to the most significant signals (the slopes are 
given too). 



smaller the damping, the more the mode occupation is 
driven out of equilibrium, i.e., to a larger average excita- 
tion level. In the extreme case of no damping 7damp = 
(dotted curve) [the externally undamped limit in Ref. F78| . 
the local heating is maximal. On the other hand, a suf- 
ficiently large damping may effectively prevent phonon 
heating [the externally damped limit in Ref. [tII . From 
Fig. [15] we see that with a phonon damping as large as 
200 fieV/h the slope has vanished. 

Figure [T3] also compares our theoretical results to the 
original experimental measurements by Agrai't et al^ 
(noisy curves). The four experimental characteristics 
(aligned with the calculated zero-bias conductance) cor- 
responds to a presumably 7-atom long gold wire under 
different states of strain recorded at low temperatures 
T = 4.2 K. From this plot it is clear that theory and ex- 
periment are in excellent agreement with respect to the 
position of the phonon signal and the magnitude of the 
dominant drop. One also notices the indication of a sec- 
ondary phonon feature below 5 meV in all curves. But 
what is particularly interesting is that the measured con- 
ductance slopes beyond the threshold seem to agree well 



with a phonon damping of the order 5-50 ^eV, which is 
further quite reasonable according to our estimate above. 
The only feature which is not perfectly reproduced is 
the experimental width of phonon signal lincshape — as 
seen from the derivative of the conductance dG/dV in 
the lower part of the figure — which is somewhat wider 
than the calculated ones (which for comparison also in- 
cludes the instrumental lock-in broadening corresponding 
Kms = 1 meV). 



V. HYDROCARBON MOLECULES BETWEEN 
GOLD CONTACTS 

The general method described in Sec. Ulland lllll is ap- 
plicable to many other systems than atomic gold wires. 
Examples of systems where it is interesting to apply this 
method include wires and contacts of other metals as 
well as individual molecules. In fact, we have already 
used the present method to study conjugated and sat- 
urated hydrocarbon molecules in between gold surfaces, 
see Ref. rza The purpose of this section is to illustrate 
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FIG. 14: (Color online) ABL-mode broadening due to cou- 
pling to bulk phonons. The spectrum Bx{e) corresponds to 
the important ABL-mode for a 7-atom wire (L = 29.20A). 
By fitting the calculated points with a Lorentzian we extract 
a full- width half max (FWHM) broadening of 27damp = S^eV 
and a frequency shift of Sujx = — 6/ieV. The inset shows the 
calculated total density of states for bulk An (full line) , as well 
as a decomposition in the direction of the transport (dashed 
red curve) and in the transverse direction (dotted blue curve). 
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that our method is general enough to apply to many sys- 
tems; especially that the LOE approximation is likely to 
be valid for a range of systems where, at first glance, it 
is not expected to work. 

We start with a brief description of our previous 
results^^ motivated by the recent experiments by Kushm- 
erick et ali^ They measured the inelastic scattering sig- 
nal through three different molecules (Cll, OPV, and 
OPE) connected to gold electrodes by means of a cryo- 
genic crossed-wire tunnel junction setup. Since the num- 
ber of molecules present in the experimentally realized 
junctions is unknown it is advantageous to look at the 
inelastic electron tunneling spectroscopy (lETS) signal 
defined as 



di/dv ■ 



(57) 



which — if the current / simply scales with the number of 
molecules — is independent of the number of molecules in 
the junction. 

In Ref. [tI, we used the present LOE method to model 
the lETS spectra for each of these three molecules. As 
an example, Fig. [1^] shows the calculated and mea- 
sured lETS spectrum in the case of the conjugated OPE 
molecule [inset of Fig. [TWb)]. It is seen that our theory 
reproduces the positions and relative heights of the in- 
elastic scattering peaks. The three main peaks are given 



FIG. 15: (Col or on line) Comparison between theory and ex- 
periment (Ref. IIITI ) for the inelastic conductance of an atomic 
gold wire. The measured characteristics (noisy black curves) 
correspond to different states of strain of wire (around 7 atoms 
long). The calculated results (smooth colored lines) are for 
the 7-atom wire at L = 29.20 A using different values for the 
external damping as indicated in the upper right corner of 
the plot (in units of fj,eV/h). The dashed curve is the cal- 
culated result in the externally undamped limit (7damp ~ 0). 
The lower plot is the numerical derivative of the conductance, 
where the experimental curves have been offset by multiples of 
Go/V for clarity. Note the indication of a secondary phonon 
feature below 5 meV in all curves. The temperature is T = 4.2 
K and the lock-in modulation voltage Kms = 1 meV (in both 
theory and experiment). 



by four types of vibrations; one type is affecting the C- 
S stretch whereas the other three involve the distortion 
of the C backbone of the molecule. In our calculation 
the region of dynamic atoms includes 54 atoms corre- 
sponding to 162 vibrational modes (18 Au surface atoms 
and 36 atoms in the molecule). We thus see that the 
lETS spectrum must be related to certain selection rules 
that describe why only a few vibrational modes affect the 
current. These selection rules may be understood from 
studying the electron scattering states and the symmetry 
of the e-ph interaction^^ For the other two molecules 
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FIG. 16: (Color online) Calculated lETS spectrum for an 
OPE molecule compared to the experimental data from 
Ref. [l^. Each of the three inelastic scattering peaks arise 
from different kinds of vibrations localized on the molecule. 
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FIG. 17: (Color online) Calculated lETS spectra for (a) an 
OPV molecule and (b) an OPE molecule. The chemical struc- 
ture of these hydrocarbon molecules are shown in the insets. 
The two plots show that the simple LOE scheme predicts the 
same lETS spectrum as the full SCBA (if one neglects the 
elastic variation). 



(OPV and Cll) wc found a similar good agreement 
with the experiments by Kushmerick et al. However, the 
transmission T{e) through these three molecules is actu- 



ally varying significantly with energy, since the electron 
conduction process involves states around the Fermi en- 
ergy that lie in the gap between the molecular levels. 
For instance, in an energy window of 0.4 eV this varia- 
tion is of the order T(£f - 0.2eV)/T(£F + 0.2eV) « 4 for 
the OPE molecule. Accordingly the use of the LOE ap- 
proximation might seem inappropriate for these systems. 
With a detailed comparison between LOE and full SCBA 
calculations (including this energy dependence) we can 
nevertheless show that the LOE approximation provides 
effectively the same prediction for the lETS spectrum. 
This comparison is found in Fig. 1171 

Since the SCBA is computationally expensive it is not 
realistic to use the same high accuracy as for LOE calcu- 
lations. We therefore reduce the device subspace and the 
region of dynamic atoms to include only the molecule. 
Furthermore we use a smaller SZP basis set describing 
the OPE (OPV) molecule reducing the device subspace 
to 264 (280) atomic orbitals. Finally we include only the 
5 (3) most important vibrational modes (selected from a 
LOE calculation). With these simplifications we calcu- 
lated the current for 81 (61) bias points using an average 
of 9 (8) iterations to converge the SCBA on an energy 
grid of approximately 500 points. These SCBA calcula- 
tions required 40 (18) hours on 10 Pentium-4 processors 
working in parallel. For comparison, the corresponding 
LOE calculations can be performed in less than 1 minute 
on a single Pentium-4 processor. 

The results shown in Fig. [T7] reveal that the LOE ap- 
proximation captures the inelastic scattering signal with 
a very satisfactory accuracy. The main discrepancy be- 
tween LOE and SCBA is directly related to the elastic 
part of the transport which can easily be corrected for 
without solving the full SCBA equations, cf. Sec. IIVDI 
We have thus used our implementation of SCBA to jus- 
tify that the simpler LOE scheme can actually be applied 
for the lETS spectra of the hydrocarbon molecules. This 
is not a trivial result because the energy variation in the 
transmission around the Fermi energy for these systems 
seems to violate one of the fundamental assumptions of 
the LOE. 



VI. CONCLUSIONS AND OUTLOOK 

In this paper we have presented a first-principles 
method for calculating the effects of vibrations and e- 
ph couplings in the electronic transport properties of an 
atomic-scale device. Our implementation that extends 
the Siesta implementation of Kohn-Sham DFT and the 
Transiesta scheme for elastic transport is described in 
detail, highlighting the important computational steps 
for the complete analysis. The inelastic transport prob- 
lem is addressed using the NEGF formalism with the 
e-ph interaction treated up to the level of SCBA. We 
also describe the computationally simple LOE scheme. 
As illustrations of the methodology we have applied it to 
model the phonon signals in the conductance of atomic 
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gold wires and hydrocarbon molecules between gold sur- 
faces. In both cases the comparison with experimental 
results is very satisfactory. While we expect our method 
to be successful for a wide range of nanoscale systems, 
there are also some important aspects where further re- 
search and development may lead to improvements. We 
therefore close this paper with an outlook of some of the 
challenges we believe are important. 

While we have argued that the vibrations for the sys- 
tems considered here are reasonably well described by 
free phonon Green's functions, there might also be sit- 
uations where the phonon system has to be treated be- 
yond free dynamics, e.g., by including self-energies from 
e-h pair damping, anharmonic phonon-phonon couplings 
(inside the device), and resonant phonon-phonon cou- 
plings (between device and electrodes). And as we have 
shown in this work, these precise damping conditions of 
the phonons are governing the device heating. Another 
issue is the bias-induced changes in geometry and e-ph 
couplings. Further development along these lines might 
thus lead to a better understanding of transport in the 
high-bias regime. On the more technical side, it would 
be interesting to extend the present scheme to describe 
the interplay between c-ph couplings and other delicate 
effects such as spin-polarized currents, spin-orbit cou- 
plings, etc. For instance, phonon heating could mediate 
an important effective interaction between the two spin 
channels. 

In conclusion, the present paper contributes to the 
evolving understanding of phonon scattering and local 
heating in nanoscale systems. These effects are impor- 
tant to elucidate the structural properties from the elec- 
tronic transport characteristics and ultimately for the 
stability of devices. 



APPENDIX A: HILBERT TRANSFORM 

The purpose of this appendix is to discuss efficient nu- 
merical ways to approximate the Hilbert transform of a 
continuous function f{x), here defined asi^ 
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(Al) 



where V denotes the Cauchy principal value integral. 

We approximate the function f{x) by a linear interpo- 
lation fi{x) to the values fi = f{xi) known at the dis- 
crete grid points {xi}. This we can write in the following 
way 



fix) - fi{x) 



N 

E 

i=l 



(A2) 



where the kernel function associated with the linear in- 
terpolation is 



■>pi{x) 



X 



Xi Xi- 



Xi+1 



-[9{xi - x) - 9{xi-i - x)] 

X 

— [9{xi+i - x) - 9{xi - x)]. 



(A3) 



On this form we implicitly assume that the function falls 
off to zero at the ends of the grid, i.e., that the function 
has finite support. We can then approximate the Hilbert 
transform of f{x) by the Hilbert transform of fi{x), i.e.. 



= -V 

TT 

N 



dx 



(A4) 



where we have identified a transformation kernel 
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dx 



ip,{x) 



Xi-1 



Xi+l 



Xi+1 - Xi 



In 



■In 



Xi-1 - 
^i+l 



(A5) 



Having determined the matrix Kji corresponding to a 
given grid {xi}, the Hilbert transform amounts to a 
matrix-vector product operation. With N grid points 
this scales as 0{N^). 

A typical situation is that of an equidistant grid Xi — 
Xi-1 = A (for all i), where a more effective algorithm can 
be devised. In this case we can write Xi — Xj = {i — j) A, 
and the kernel function, that becomes a function of the 
index difference m — j i only, reduces to 



— \ — (m — 1) ln(m — 1) 
+2m Inm — (m + 1) ln(m 



1)]. (A6) 
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1 KfL Ji has 



The Hilbert transform Ti-xifiyixj) = 
then taken the form of a discrete convolution which effec- 
tively can be calculated with the Fast Fourier transform 
(FFT) algorithm. This scales only as O(A^lniV). 



APPENDIX B: NUMERICAL 
IMPLEMENTATION OF SCBA 

Calculating the current numerically using the SCBA 
is highly nontrivial for large systems. This appendix dis- 
cusses our solutions to the main difficulties encountered 
within the SCBA. We exemplify the size and scope of 
the calculations, e.g., the sizes of matrices and the en- 
ergy grid, with values taken from the SCBA calculation 
presented in Sec. |V|on the OPE molecule. 

The current and power expressions Eq. (PT]) and ([5^ 
arc integrated numerically using a third order polynomial 
interpolation. Since the inelastic signal is typically small, 
the current has to be determined with a high accuracy, 
which implies a fine resolution of the energy grid for the 
integration. Further, the range of this grid has to include 
not only the bias window but also additional energies due 
to the nonlocal character (in energy) of the Hilbert trans- 
form, cf. Eq. (|39p . These limitations make a nonuniform 
grid preferable. We thus construct a dense grid around 
each of the important energies e = ^l,r, ^J■L,R ± fi^A, • ■ • 
and a coarser one elsewhere. The resolution of the fine 
grid is determined by the temperature and should have 
a point separation around Se < 0.5 keT. For the OPE 
molecule we found it adequate at T = 40 K to use a fine 
grid with 5e = 1.7 meV and a coarse grid with Ae = 10.0 
mcV spanning the energy range [—0.5,0.5] cV. With a 
nonuniform grid the necessary number of energy points 
may thus be reduced. 

The solution of the SCBA approximation requires sub- 
stantial amounts of CPU time and memory. Analyzing 
the memory requirements we find that we need to retain 
G5''-(e) and Sj^h^e) 

in memory. Each of these matrices 

requires a memory allocation of ©(A'giid A^^asis) bytes, 
where A^grid is the number of grid points, and A^basis 
the size of the electronic basis. For the OPE calcula- 
tion in Sec. |V] each matrix takes up 500 Megabytes of 
memory (500 energy points x 250^ matrix size x 16 
bytes/complex number). In addition to the demand- 
ing memory requirement, significant computational time 
(400 CPU hours in total) is needed. 

The computationally heaviest part is the calculation of 
Eq. pS]) . which we rewrite as 

Sg^(£) = ^M^[(nA)G^(£±nc.A) (Bl) 

A 

From this equation we see that the CPU time scales 
as ^^(Afph A^grid Af basis -^itcr) [siucc cach matrix multipli- 
cation scales as £>( A^basis)]; where A^ph is the number 



of vibrational modes and Alitor the number of iterations 
needed for self-consistency of the SCBA. 

We have overcome the memory and computational re- 
quirements by a parallelization of our computer code by 
dividing the energy grid over the available processors. 
The only significant complication is the evaluation of 
Eq. (|Bip , where quantities couple across the energy divi- 
sion. To overcome this, wc first redistribute the Green's 
functions G^(e) over the processors by changing from en- 
ergy division to matrix indices division. Then the energy- 
shifted Green's functions can be added for each matrix 
index. Next we transform the outcome back to energy 
division and carry out the matrix multiplications with 
Ma- We have implemented this procedure efficiently in 
a way that lets the necessary communication occur while 
other calculations are running, i.e., while the lesser part 
of the equation is being communicated between proces- 
sors, the matrix multiplications for the greater part are 
being computed and vice versa. In practice, this par- 
allelization works very well and the computation time 
scales almost linearly with the number of processors. 



APPENDIX C: SIGNAL BROADENING BY 
LOCK-IN MODULATION VOLTAGE 



As discussed in Sec. IIIIFI the lock-in technique for 
measuring the differential conductance (and derivatives) 
introduces a broadening of the intrinsic current-voltage 
characteristics due to a finite modulation voltage. The 
basic idea is to measure the frequency components of the 
current at multiples of the applied harmonic modulation, 
since these relates to the derivatives of the current. Fol- 
lowing Hansmaj^^ we can analytically write the frequency 
components as the following averages over an oscillation 
period 



2tt/u 



I,., = — J I[V + Acos{ujt)]cos{ujt)dt 





2 f dI{V + Ax) 
dV 



\/l^-x^ dx, 



(CI) 



and 



I[V + Acosiujt)] cos{2ujt) dt 



8 f d^I{V + Ax) , 2N3/2 



dV^ 



\-xy dx, (C2) 



where the modulation amplitude is A = v^Umis- The 
partial integrations carried out above show that the com- 
ponents I^^ and l2ui are convolutions of the exact first and 
second derivatives of the current with certain functions 
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proportional to \J\ — and (1 — a;^)'^/^, respectively. If 
we assume that the inelastic signal has no intrinsic width, 
the inelastic conductance change is proportional to a step 
function d{eV—hujx) and the second derivative to a delta 



function S{eV—HLLj\). With these functional forms the in- 
tegrals can be evaluated, leading to a modulation broad- 
ening of the first (second) derivative of approximately 

2.45 Vrn.s (1.72 V,n,s)- 
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